A Model for Understanding Numerical Stability

We present a model of roundoff error analysis that combines simplicity with predictive power. Though not considering all sources of roundoff within an algorithm, the model is related to a recursive roundoff error analysis and therefore capable of correctly predicting stability or instability of an algorithm. By means of nontrivial examples, such as the componentwise backward stability analysis of Gaussian elimination with a single iterative refinement step, we demonstrate that the model even yields quantitative backward error bounds that show all the known problem-dependent terms (with the exception of dimension-dependent constants, which are the weak spot of any a priori analysis). The model can serve as a convenient tool for teaching or as a heuristic device to discover stability results before entering a further, detailed analysis

On the Scaling Limits of Determinantal Point Processes with Kernels Induced by Sturm-Liouville Operators

By applying an idea of Borodin and Olshanski [J. Algebra 313 (2007), 40-60], we study various scaling limits of determinantal point processes with trace class projection kernels given by spectral projections of selfadjoint Sturm-Liouville operators. Instead of studying the convergence of the kernels as functions, the method directly addresses the strong convergence of the induced integral operators. We show that, for this notion of convergence, the Dyson, Airy, and Bessel kernels are universal in the bulk, soft-edge, and hard-edge scaling limits. This result allows us to give a short and unified derivation of the known formulae for the scaling limits of the classical random matrix ensembles with unitary invariance, that is, the Gaussian unitary ensemble (GUE), the Wishart or Laguerre unitary ensemble (LUE), and the MANOVA (multivariate analysis of variance) or Jacobi unitary ensemble (JUE)

On the convergence rates of Gauss and Clenshaw-Curtis quadrature for functions of limited regularity

We study the optimal general rate of convergence of the n-point quadrature rules of Gauss and Clenshaw-Curtis when applied to functions of limited regularity: if the Chebyshev coefficients decay at a rate O(n^{-s-1}) for some s > 0, Clenshaw-Curtis and Gauss quadrature inherit exactly this rate. The proof (for Gauss, if 0 < s < 2, there is numerical evidence only) is based on work of Curtis, Johnson, Riess, and Rabinowitz from the early 1970s and on a refined estimate for Gauss quadrature applied to Chebyshev polynomials due to Petras (1995). The convergence rate of both quadrature rules is up to one power of n better than polynomial best approximation; hence, the classical proof strategy that bounds the error of a quadrature rule with positive weights by polynomial best approximation is doomed to fail in establishing the optimal rate.Comment: 7 pages, the figure of the revision has an unsymmetric example, to appear in SIAM J. Numer. Ana

Sparsing in Real Time Simulation

Modelling of mechatronical systems often leads to large DAEs with stiff components. In real time simulation neither implicit nor explicit methods can cope with such systems in an efficient way: explicit methods have to employ too small steps and implicit methods have to solve too large systems of equations. A solution of this general problem is to use a method that allows manipulations of the Jacobian by computing only those parts that are necessary for the stability of the method. Specifically, manipulation by sparsing aims at zeroing out certain elements of the Jacobian leading t a structure that can be exploited using sparse matrix techniques. The elements to be neglected are chosen by an a priori analysis phase that can be accomplished before the real-time simulaton starts. In this article a sparsing criterion for the linearly implicit Euler method is derived that is based on block diagnonalization and matrix perturbation theory

Automatic Deformation of Riemann-Hilbert Problems with Applications to the Painlev\'e II Transcendents

The stability and convergence rate of Olver's collocation method for the numerical solution of Riemann-Hilbert problems (RHPs) is known to depend very sensitively on the particular choice of contours used as data of the RHP. By manually performing contour deformations that proved to be successful in the asymptotic analysis of RHPs, such as the method of nonlinear steepest descent, the numerical method can basically be preconditioned, making it asymptotically stable. In this paper, however, we will show that most of these preconditioning deformations, including lensing, can be addressed in an automatic, completely algorithmic fashion that would turn the numerical method into a black-box solver. To this end, the preconditioning of RHPs is recast as a discrete, graph-based optimization problem: the deformed contours are obtained as a system of shortest paths within a planar graph weighted by the relative strength of the jump matrices. The algorithm is illustrated for the RHP representing the Painlev\'e II transcendents.Comment: 20 pages, 16 figure

The Singular Values of the GOE

As a unifying framework for examining several properties that nominally involve eigenvalues, we present a particular structure of the singular values of the Gaussian orthogonal ensemble (GOE): the even-location singular values are distributed as the positive eigenvalues of a Gaussian ensemble with chiral unitary symmetry (anti-GUE), while the odd-location singular values, conditioned on the even-location ones, can be algebraically transformed into a set of independent $\chi$-distributed random variables. We discuss three applications of this structure: first, there is a pair of bidiagonal square matrices, whose singular values are jointly distributed as the even- and odd-location ones of the GOE; second, the magnitude of the determinant of the GOE is distributed as a product of simple independent random variables; third, on symmetric intervals, the gap probabilities of the GOE can be expressed in terms of the Laguerre unitary ensemble (LUE). We work specifically with matrices of finite order, but by passing to a large matrix limit, we also obtain new insight into asymptotic properties such as the central limit theorem of the determinant or the gap probabilities in the bulk-scaling limit. The analysis in this paper avoids much of the technical machinery (e.g. Pfaffians, skew-orthogonal polynomials, martingales, Meijer $G$-function, etc.) that was previously used to analyze some of the applications.Comment: Introduction extended, typos corrected, reference added. 31 pages, 1 figur

Asymptotic expansions relating to the lengths of longest monotone subsequences of involutions

We study the distribution of the length of longest monotone subsequences in random (fixed-point free) involutions of $n$ integers as $n$ grows large, establishing asymptotic expansions in powers of $n^{-1/6}$ in the general case and in powers of $n^{-1/3}$ in the fixed-point free cases. Whilst the limit laws were shown by Baik and Rains to be one of the Tracy-Widom distributions $F_\beta$ for $\beta=1$ or $\beta=4$, we find explicit analytic expressions of the first few finite-size correction terms as linear combinations of higher order derivatives of $F_\beta$ with rational polynomial coefficients. Our derivation is based on a concept of generalized analytic de-Poissonization and is subject to the validity of certain hypotheses for which we provide compelling (computational) evidence. In a preparatory step expansions of the hard-to-soft edge transition laws of L$\beta$E are studied, which are lifted into expansions of the generalized Poissonized length distributions for large intensities. (This paper continues our work arXiv:2301.02022, which established similar results in the case of general permutations and $\beta=2$.)Comment: 44 pages, 5 figure, 3 table

A Jentzsch-Theorem for Kapteyn, Neumann, and General Dirichlet Series

Comparing phase plots of truncated series solutions of Kepler's equation by Lagrange's power series with those by Bessel's Kapteyn series strongly suggest that a Jentzsch-type theorem holds true not only for the former but also for the latter series: each point of the boundary of the domain of convergence in the complex plane is a cluster point of zeros of sections of the series. We prove this result by studying properties of the growth function of a sequence of entire functions. For series, this growth function is computable in terms of the convergence abscissa of an associated general Dirichlet series. The proof then extends, besides including Jentzsch's classical result for power series, to general Dirichlet series, to Kapteyn, and to Neumann series of Bessel functions. Moreover, sections of Kapteyn and Neumann series generally exhibit zeros close to the real axis which can be explained, including their asymptotic linear density, by the theory of the distribution of zeros of entire functions.Comment: corrected Hankel asymptotics on p. 14; 16 pages, 4 figure